IMPORTY BIBLIOTEK¶

In [1]:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import geopandas as gpd
import seaborn as sns
import random
from statsmodels.graphics.tsaplots import plot_pacf

ZAPOZNANIE Z DANYMI¶

In [2]:
data = pd.read_excel(r"C:\Users\User\Desktop\magister\semestr1\ML\projekt_ocena_2\2019_PM25_1g.xlsx", skiprows = 5)
data.rename(columns = {'Kod stanowiska':'Date'}, inplace = True) 
# data.set_index('Date', inplace = True)
# data.index = pd.to_datetime(data.index)
first_column = data.iloc[:, 0]  
data = data.drop(data.columns[0], axis = 1) 
data.columns = data.columns.str.replace('-PM2.5-1g', '')
data.insert(0, 'Kod stanowiska', first_column)
data.rename(columns = {'Kod stanowiska': 'Date'}, inplace = True)
data.head(5)
C:\Users\User\AppData\Local\Temp\ipykernel_23260\4022153991.py:7: FutureWarning: The default value of regex will change from True to False in a future version.
  data.columns = data.columns.str.replace('-PM2.5-1g', '')
Out[2]:
Date DsDusznikMOB DsJaworMOB DsJelGorOgin DsWrocAlWisn DsWrocWybCon KpBydPlPozna KpMogiNowMOB KpToruDziewu KpWloclOkrze ... SlBielPartyz SlKatoKossut SlZlotPotLes WmElbBazynsk WmGoldUzdrowMOB WmOlsPuszkin WpKaliSawick ZpSzczAndr01 ZpSzczBudzWosMOB ZpSzczPils02
0 2019-01-01 01:00:00 33.40530 51.38780 118.7730 102.0900 107.0610 64.1177 NaN 24.0030 51.317 ... 110.1990 76.7306 26.3444 34.3706 14.9449 40.9183 75.2000 NaN NaN 73.93500
1 2019-01-01 02:00:00 13.80280 28.49950 110.0640 63.6111 55.9187 43.8401 NaN 33.6542 30.698 ... 73.4132 54.4664 19.0619 23.1494 10.7420 25.9358 47.9076 NaN NaN 11.78830
2 2019-01-01 03:00:00 9.94056 11.12060 107.9410 48.3540 41.3488 22.8383 NaN 13.6030 28.262 ... 50.2355 50.4599 43.7717 21.0711 12.0391 24.5725 22.8309 5.57095 NaN 8.69917
3 2019-01-01 04:00:00 6.75889 5.57358 94.5489 34.6621 29.8697 20.1829 NaN 17.4302 26.522 ... 37.5872 34.8090 64.0139 21.1671 13.1849 20.6336 20.5900 5.77369 NaN 5.96861
4 2019-01-01 05:00:00 7.88722 6.56224 67.8800 14.2870 17.6000 18.7345 NaN 23.0878 24.260 ... 22.6446 30.6517 43.6111 21.0774 14.0005 19.4194 27.0838 6.15494 NaN 7.80778

5 rows × 64 columns

Zdecydowano się znaleźć dane dotyczące lokalizacji powyższych stacji pomiarowych, aby dodać uwzględnienie czynnika przestrzennego w dalszej analizie. Być może będą przydatne m.in. przy uzupełnianiu brakujących wartości lub przy sprawdzaniu zależności pomiędzy pomiarami w danym województwie. Zmieniono nazwy kolumn na takie, które odpowiadają nazwom stacji.

In [3]:
kody_stacji = pd.read_excel(r"C:\Users\User\Desktop\magister\semestr1\ML\projekt_ocena_2\Metadane oraz kody stacji i stanowisk pomiarowych (1).xlsx")
kody_stacji.head(5)
Out[3]:
Nr Kod stacji Kod międzynarodowy Nazwa stacji Stary Kod stacji \n(o ile inny od aktualnego) Data uruchomienia Data zamknięcia Typ stacji Typ obszaru Rodzaj stacji Województwo Miejscowość Adres WGS84 φ N WGS84 λ E
0 1 DsBialka NaN Białka NaN 1990-01-03 2005-12-31 przemysłowa podmiejski kontenerowa stacjonarna DOLNOŚLĄSKIE Białka NaN 51.197783 16.117390
1 2 DsBielGrot NaN Bielawa - ul. Grota Roweckiego NaN 1994-01-02 2003-12-31 tło miejski w budynku DOLNOŚLĄSKIE Bielawa ul. Grota Roweckiego 6 50.682510 16.617348
2 3 DsBogatFrancMOB PL0602A Bogatynia Mobil DsBogatMob 2015-01-01 2015-12-31 tło miejski mobilna DOLNOŚLĄSKIE Bogatynia ul. Francuska/Kręta 50.940998 14.916790
3 4 DsBogChop PL0315A Bogatynia - Chopina NaN 1996-01-01 2013-12-31 przemysłowa miejski kontenerowa stacjonarna DOLNOŚLĄSKIE Bogatynia ul. Chopina 35 50.905856 14.967175
4 5 DsBogZatonieMob PL0576A Bogatynia - Mobil NaN 2012-01-01 2012-12-31 przemysłowa miejski mobilna DOLNOŚLĄSKIE Bogatynia ul. Konrada, Zatonie 50.943245 14.913327
In [4]:
kody_stacji = kody_stacji[['Kod stacji','Nazwa stacji', 'Stary Kod stacji \n(o ile inny od aktualnego)', 'Typ obszaru', 'Województwo', 'Miejscowość', 'WGS84 φ N', 'WGS84 λ E']]
kody_stacji.rename(columns = {'WGS84 φ N': 'lat'}, inplace = True)
kody_stacji.rename(columns = {'WGS84 λ E': 'lon'}, inplace = True)
kody_stacji.head(5)
Out[4]:
Kod stacji Nazwa stacji Stary Kod stacji \n(o ile inny od aktualnego) Typ obszaru Województwo Miejscowość lat lon
0 DsBialka Białka NaN podmiejski DOLNOŚLĄSKIE Białka 51.197783 16.117390
1 DsBielGrot Bielawa - ul. Grota Roweckiego NaN miejski DOLNOŚLĄSKIE Bielawa 50.682510 16.617348
2 DsBogatFrancMOB Bogatynia Mobil DsBogatMob miejski DOLNOŚLĄSKIE Bogatynia 50.940998 14.916790
3 DsBogChop Bogatynia - Chopina NaN miejski DOLNOŚLĄSKIE Bogatynia 50.905856 14.967175
4 DsBogZatonieMob Bogatynia - Mobil NaN miejski DOLNOŚLĄSKIE Bogatynia 50.943245 14.913327

Wczytano dane znalezione na stronie GIOS, zawierające szczegółowe informacje o stacjach - ich lokalizację, województwo, typ obszaru itd. Wybrano kolumny, które uważano za istotne w dalszej analizie. Zmieniono nazwy kolumn na przyjemniejsze w użytku.

In [5]:
data_melted = data.melt(id_vars = 'Date', var_name = 'Kod stacji', value_name='PM2.5') 
data_melted.head(10)
Out[5]:
Date Kod stacji PM2.5
0 2019-01-01 01:00:00 DsDusznikMOB 33.40530
1 2019-01-01 02:00:00 DsDusznikMOB 13.80280
2 2019-01-01 03:00:00 DsDusznikMOB 9.94056
3 2019-01-01 04:00:00 DsDusznikMOB 6.75889
4 2019-01-01 05:00:00 DsDusznikMOB 7.88722
5 2019-01-01 06:00:00 DsDusznikMOB 14.40720
6 2019-01-01 07:00:00 DsDusznikMOB 7.48833
7 2019-01-01 08:00:00 DsDusznikMOB 5.67083
8 2019-01-01 09:00:00 DsDusznikMOB 11.72390
9 2019-01-01 10:00:00 DsDusznikMOB 11.72530

Zdecydowano się przekształcić DataFrame 'data', tak aby usprawnić łączenie ze szczegółowymi kodami stacji. Teraz nazwy stacji nie są w osobnych kolumnach tylko w jednej połączonej.

In [6]:
data_merged = pd.merge(data_melted, kody_stacji, on = 'Kod stacji', how = 'left')
data_merged.head(5)
Out[6]:
Date Kod stacji PM2.5 Nazwa stacji Stary Kod stacji \n(o ile inny od aktualnego) Typ obszaru Województwo Miejscowość lat lon
0 2019-01-01 01:00:00 DsDusznikMOB 33.40530 Duszniki-Zdrój NaN miejski DOLNOŚLĄSKIE Duszniki-Zdrój 50.402645 16.393319
1 2019-01-01 02:00:00 DsDusznikMOB 13.80280 Duszniki-Zdrój NaN miejski DOLNOŚLĄSKIE Duszniki-Zdrój 50.402645 16.393319
2 2019-01-01 03:00:00 DsDusznikMOB 9.94056 Duszniki-Zdrój NaN miejski DOLNOŚLĄSKIE Duszniki-Zdrój 50.402645 16.393319
3 2019-01-01 04:00:00 DsDusznikMOB 6.75889 Duszniki-Zdrój NaN miejski DOLNOŚLĄSKIE Duszniki-Zdrój 50.402645 16.393319
4 2019-01-01 05:00:00 DsDusznikMOB 7.88722 Duszniki-Zdrój NaN miejski DOLNOŚLĄSKIE Duszniki-Zdrój 50.402645 16.393319

Połączono dwa zestawy w jeden szczegółowy.

In [7]:
not_merged = data_merged[data_merged['Nazwa stacji'].isnull()]
not_merged['Kod stacji'].unique()
Out[7]:
array(['LbNaleczow', 'MzKonJezMos', 'PdSuwPulaskp', 'PmGdaLeczk08',
       'ZpSzczAndr01', 'ZpSzczPils02'], dtype=object)

6 nazw stacji nie znalzało swoich odpowiedników, więc zostaną podjęte próby ich odnalezienia, o ile istnieją

In [8]:
kody_stacji[kody_stacji['Stary Kod stacji \n(o ile inny od aktualnego)'] == 'LbNaleczow']
Out[8]:
Kod stacji Nazwa stacji Stary Kod stacji \n(o ile inny od aktualnego) Typ obszaru Województwo Miejscowość lat lon
In [9]:
kody_stacji_notna = kody_stacji[kody_stacji['Stary Kod stacji \n(o ile inny od aktualnego)'].notna()]
kody_stacji_notna[kody_stacji_notna['Stary Kod stacji \n(o ile inny od aktualnego)'].str.contains('LbNaleczow')]
Out[9]:
Kod stacji Nazwa stacji Stary Kod stacji \n(o ile inny od aktualnego) Typ obszaru Województwo Miejscowość lat lon
296 LbNaleczAlMa Nałęczów, al. Małachowskiego LbNaleczowMOB, LbNaleczow podmiejski LUBELSKIE Nałęczów 51.284931 22.210242
In [10]:
kody_stacji[kody_stacji['Stary Kod stacji \n(o ile inny od aktualnego)'] == 'MzKonJezMos']
Out[10]:
Kod stacji Nazwa stacji Stary Kod stacji \n(o ile inny od aktualnego) Typ obszaru Województwo Miejscowość lat lon
536 MzKonJezWieMOB Konstancin-Jeziorna, ul. Wierzejewskiego MzKonJezMos podmiejski MAZOWIECKIE Konstancin-Jeziorna 52.080625 21.111186
In [11]:
kody_stacji[kody_stacji['Stary Kod stacji \n(o ile inny od aktualnego)'] == 'PdSuwPulaskp']
Out[11]:
Kod stacji Nazwa stacji Stary Kod stacji \n(o ile inny od aktualnego) Typ obszaru Województwo Miejscowość lat lon
688 PdSuwPulask2 Suwałki, ul. Pułaskiego 26 PdSuwPulaskp miejski PODLASKIE Suwałki 54.115897 22.938464
In [12]:
kody_stacji[kody_stacji['Stary Kod stacji \n(o ile inny od aktualnego)'] == 'PmGdaLeczk08']
Out[12]:
Kod stacji Nazwa stacji Stary Kod stacji \n(o ile inny od aktualnego) Typ obszaru Województwo Miejscowość lat lon
In [13]:
kody_stacji_notna[kody_stacji_notna['Stary Kod stacji \n(o ile inny od aktualnego)'].str.contains('PmGdaLeczk08')]
Out[13]:
Kod stacji Nazwa stacji Stary Kod stacji \n(o ile inny od aktualnego) Typ obszaru Województwo Miejscowość lat lon
799 PmGdaLeczkow Gdańsk, ul. Leczkowa Pm.a08a, PmGdaLeczk08, PmGdaLecz08m miejski POMORSKIE Gdańsk 54.380279 18.620274
In [14]:
kody_stacji[kody_stacji['Stary Kod stacji \n(o ile inny od aktualnego)'] == 'ZpSzczAndr01']
Out[14]:
Kod stacji Nazwa stacji Stary Kod stacji \n(o ile inny od aktualnego) Typ obszaru Województwo Miejscowość lat lon
In [15]:
kody_stacji_notna[kody_stacji_notna['Stary Kod stacji \n(o ile inny od aktualnego)'].str.contains('ZpSzczAndr01')]
Out[15]:
Kod stacji Nazwa stacji Stary Kod stacji \n(o ile inny od aktualnego) Typ obszaru Województwo Miejscowość lat lon
1082 ZpSzczAndrze Szczecin, ul. Andrzejewskiego ZpSzczecin001, ZpSzczAndr01 miejski ZACHODNIOPOMORSKIE Szczecin 53.380975 14.663347
In [16]:
kody_stacji[kody_stacji['Stary Kod stacji \n(o ile inny od aktualnego)'] == 'ZpSzczPils02']
Out[16]:
Kod stacji Nazwa stacji Stary Kod stacji \n(o ile inny od aktualnego) Typ obszaru Województwo Miejscowość lat lon
In [17]:
kody_stacji_notna[kody_stacji_notna['Stary Kod stacji \n(o ile inny od aktualnego)'].str.contains('ZpSzczPils02')]
Out[17]:
Kod stacji Nazwa stacji Stary Kod stacji \n(o ile inny od aktualnego) Typ obszaru Województwo Miejscowość lat lon
1090 ZpSzczPilsud Szczecin, ul. Piłsudskiego ZpSzczecin002, ZpSzczPils02 miejski ZACHODNIOPOMORSKIE Szczecin 53.432169 14.5539
In [18]:
data_melted['Kod stacji'] = data_merged['Kod stacji'].replace({'LbNaleczow' : 'LbNaleczAlMa',
                                                                 'MzKonJezMos': 'MzKonJezWieMOB',
                                                                 'PdSuwPulaskp': 'PdSuwPulask2', 
                                                                 'PdSuwPulaskp': 'PdSuwPulask2', 
                                                                 'PmGdaLeczk08': 'PmGdaLeczkow',
                                                                 'ZpSzczAndr01': 'ZpSzczAndrze', 
                                                                 'ZpSzczPils02': 'ZpSzczPilsud'})

Zatąpiono stare kody stacji nowymi, aby wszystkie dane się odpowiednio połączyły.

In [19]:
data_merged = pd.merge(data_melted, kody_stacji, on = 'Kod stacji', how = 'left')
data_merged = data_merged[['Date', 'Kod stacji', 'PM2.5', 'Nazwa stacji', 'Typ obszaru', 'Województwo', 'Miejscowość', 'lat', 'lon']]
data_merged.head(5)
Out[19]:
Date Kod stacji PM2.5 Nazwa stacji Typ obszaru Województwo Miejscowość lat lon
0 2019-01-01 01:00:00 DsDusznikMOB 33.40530 Duszniki-Zdrój miejski DOLNOŚLĄSKIE Duszniki-Zdrój 50.402645 16.393319
1 2019-01-01 02:00:00 DsDusznikMOB 13.80280 Duszniki-Zdrój miejski DOLNOŚLĄSKIE Duszniki-Zdrój 50.402645 16.393319
2 2019-01-01 03:00:00 DsDusznikMOB 9.94056 Duszniki-Zdrój miejski DOLNOŚLĄSKIE Duszniki-Zdrój 50.402645 16.393319
3 2019-01-01 04:00:00 DsDusznikMOB 6.75889 Duszniki-Zdrój miejski DOLNOŚLĄSKIE Duszniki-Zdrój 50.402645 16.393319
4 2019-01-01 05:00:00 DsDusznikMOB 7.88722 Duszniki-Zdrój miejski DOLNOŚLĄSKIE Duszniki-Zdrój 50.402645 16.393319
In [20]:
not_merged = data_merged[data_merged['Nazwa stacji'].isnull()]
not_merged
Out[20]:
Date Kod stacji PM2.5 Nazwa stacji Typ obszaru Województwo Miejscowość lat lon

Upewniono się, że wszystkie dane zostały połączone.

WIZUALIZACJA¶

In [21]:
# poland = gpd.read_file(r"C:\Users\User\Desktop\magister\semestr1\ML\projekt_ocena_2\polska (1)\polska.shp")
woj = gpd.read_file(r"C:\Users\User\Desktop\magister\semestr1\ML\projekt_ocena_2\wojewodztwa\wojewodztwa.shp")
crs_polska = 'epsg:4326'
woj = woj.to_crs(crs_polska)
fig, ax = plt.subplots(figsize = (10, 10))
woj.plot(ax = ax, facecolor = 'lightblue', edgecolor = 'black')
ax.scatter(x = data_merged['lon'], y = data_merged['lat'], color = 'red', label = 'Stacje', s = 4)
plt.show()

ZASTĘPOWANIE BRAKUJĄCYCH DANYCH¶

In [22]:
(data_merged['PM2.5'].isna().groupby(data_merged['Kod stacji']).mean() * 100).sort_values(ascending=False)[0:20]
Out[22]:
Kod stacji
PmGdaLeczkow        53.059361
ZpSzczBudzWosMOB    18.493151
KpToruDziewu        17.659817
WmElbBazynsk        12.043379
KpBydPlPozna         9.691781
SkSkarZielnaMOB      9.143836
LdLodzGdansk         7.728311
PkHorZdrParkMOB      7.682648
LuZielKrotka         7.100457
WmGoldUzdrowMOB      6.929224
LuNowaSolMOB         6.803653
OpKKozBSmial         6.735160
SlKatoKossut         6.495434
LbNaleczAlMa         6.495434
MzWarWokalna         6.289954
PkPrzemGrunw         5.878995
LuWsKaziWiel         5.673516
LdZgieMielcz         5.616438
MzWarTolstoj         5.182648
KpMogiNowMOB         5.034247
Name: PM2.5, dtype: float64

Sprawdzono % brakujących wartości stężeń PM2.5 dla każdej ze stacji. Dla jednej stacji braki wynoszą ponad 50%. Braki w kolejnych stacjach są znacznie mniejsze. W 4 stacjach większe od 10%.

In [23]:
fig, axes = plt.subplots(1, 1, figsize = (20,10))

plt.plot(data_merged['Date'].unique(), data_merged['PM2.5'][data_merged['Kod stacji'] == 'PmGdaLeczkow'])
plt.xlabel('Date')
plt.ylabel('PM2.5')
plt.xlim(data_merged['Date'].min())
plt.title('PM2.5 values for station PmGdaLeczkow')
plt.show()
In [24]:
data_merged[data_merged['Województwo'] == 'PmGdaLeczkow']
Out[24]:
Date Kod stacji PM2.5 Nazwa stacji Typ obszaru Województwo Miejscowość lat lon

Z racji braku innych obserwacji w województwie pomorskim i dużej ilości brakujących danych w początkowym okresie pomiarów zdecydowano się usunąć stację.

In [25]:
data_merged = data_merged[data_merged['Kod stacji'] != 'PmGdaLeczkow']
In [26]:
data_merged.groupby('Województwo')['Kod stacji'].nunique().sort_values()
Out[26]:
Województwo
WIELKOPOLSKIE           1
LUBELSKIE               2
OPOLSKIE                2
PODLASKIE               3
WARMIŃSKO-MAZURSKIE     3
ZACHODNIOPOMORSKIE      3
ŁÓDZKIE                 3
ŚLĄSKIE                 3
ŚWIĘTOKRZYSKIE          3
KUJAWSKO-POMORSKIE      4
LUBUSKIE                4
MAŁOPOLSKIE             4
DOLNOŚLĄSKIE            5
PODKARPACKIE            7
MAZOWIECKIE            15
Name: Kod stacji, dtype: int64
In [27]:
srednia_wojewodzka_godz = data_merged.groupby(['Date', 'Województwo'])['PM2.5'].transform('mean')
data_merged['PM2.5'] = data_merged['PM2.5'].fillna(srednia_wojewodzka_godz)
data_merged['srednia_wojewodzka_godz'] = srednia_wojewodzka_godz
data_merged['PM2.5'] = data_merged['PM2.5'].fillna(srednia_wojewodzka_godz)
In [28]:
data_merged.head(5)
Out[28]:
Date Kod stacji PM2.5 Nazwa stacji Typ obszaru Województwo Miejscowość lat lon srednia_wojewodzka_godz
0 2019-01-01 01:00:00 DsDusznikMOB 33.40530 Duszniki-Zdrój miejski DOLNOŚLĄSKIE Duszniki-Zdrój 50.402645 16.393319 82.543420
1 2019-01-01 02:00:00 DsDusznikMOB 13.80280 Duszniki-Zdrój miejski DOLNOŚLĄSKIE Duszniki-Zdrój 50.402645 16.393319 54.379220
2 2019-01-01 03:00:00 DsDusznikMOB 9.94056 Duszniki-Zdrój miejski DOLNOŚLĄSKIE Duszniki-Zdrój 50.402645 16.393319 43.740992
3 2019-01-01 04:00:00 DsDusznikMOB 6.75889 Duszniki-Zdrój miejski DOLNOŚLĄSKIE Duszniki-Zdrój 50.402645 16.393319 34.282634
4 2019-01-01 05:00:00 DsDusznikMOB 7.88722 Duszniki-Zdrój miejski DOLNOŚLĄSKIE Duszniki-Zdrój 50.402645 16.393319 22.843292

Zdecydowano się zastąpić braki w danych średnią wojewódzką godzinową dla każdego z województw. Dane wciąż zawierają brakujące wartości ze względu na brak danych dla niektórych godzin.

In [29]:
(data_merged['PM2.5'].isna().groupby(data_merged['Kod stacji']).mean() * 100).sort_values(ascending=False)[0:20]
Out[29]:
Kod stacji
WpKaliSawick    2.420091
LdLodzGdansk    0.034247
LbLubObywate    0.034247
LbNaleczAlMa    0.034247
LdLodzCzerni    0.034247
LdZgieMielcz    0.034247
MzWarTolstoj    0.011416
MzPlocMiReja    0.011416
MzRadTochter    0.011416
MzSiedKonars    0.011416
MzWarAlNiepo    0.011416
SkKielTargow    0.011416
MzWarChrosci    0.011416
MzWarKondrat    0.011416
MzWarBajkowa    0.011416
MzWarWokalna    0.011416
MzOtwoBrzozo    0.011416
MzZyraRoosev    0.011416
OpKKozBSmial    0.011416
OpPrudPodgor    0.011416
Name: PM2.5, dtype: float64
In [30]:
srednia_dzienna = data_merged.groupby('Date')['PM2.5'].transform('mean')
maska_braki = (data_merged['Kod stacji'] == 'WpKaliSawick') & data_merged['PM2.5'].isna()
data_merged.loc[maska_braki, 'PM2.5']
indeksy = data_merged.loc[maska_braki, 'PM2.5'].index
# np.sum(srednia_dzienna_bez_WpKaliSawick.index == 516857)
srednia_dzienna.loc[indeksy]
data_merged.loc[maska_braki, 'PM2.5'] = srednia_dzienna.loc[indeksy]

Dla kodu stacji 'WpKaliSawick' wybrano metodę zastępowania braków średnią ze wszystkich województw.

In [31]:
(data_merged['PM2.5'].isna().groupby(data_merged['Kod stacji']).mean() * 100).sort_values(ascending=False)[0:20]
Out[31]:
Kod stacji
LdLodzGdansk       0.034247
LbLubObywate       0.034247
LbNaleczAlMa       0.034247
LdLodzCzerni       0.034247
LdZgieMielcz       0.034247
MzWarTolstoj       0.011416
MzPiasPulask       0.011416
MzPlocMiReja       0.011416
MzRadTochter       0.011416
MzSiedKonars       0.011416
MzWarAlNiepo       0.011416
SkKielTargow       0.011416
MzWarChrosci       0.011416
MzWarKondrat       0.011416
MzWarBajkowa       0.011416
MzMinMazKaziMOB    0.011416
MzWarWokalna       0.011416
MzZyraRoosev       0.011416
OpKKozBSmial       0.011416
OpPrudPodgor       0.011416
Name: PM2.5, dtype: float64
In [32]:
x = data_merged[data_merged['Kod stacji'] == 'LdZgieMielcz']['PM2.5']
x[x.isna()]
Out[32]:
117434   NaN
117435   NaN
121056   NaN
Name: PM2.5, dtype: float64
In [33]:
data_merged['PM2.5'] = data_merged.groupby('Kod stacji')['PM2.5'].apply(lambda x: x.fillna(method='ffill', limit=2).fillna(method='bfill', limit=2))
C:\Users\User\AppData\Local\Temp\ipykernel_23260\3740211283.py:1: FutureWarning: Not prepending group keys to the result index of transform-like apply. In the future, the group keys will be included in the index, regardless of whether the applied function returns a like-indexed object.
To preserve the previous behavior, use

	>>> .groupby(..., group_keys=False)

To adopt the future behavior and silence this warning, use 

	>>> .groupby(..., group_keys=True)
  data_merged['PM2.5'] = data_merged.groupby('Kod stacji')['PM2.5'].apply(lambda x: x.fillna(method='ffill', limit=2).fillna(method='bfill', limit=2))

Dla powyższych stacji braki występują tylko dla jednej, bądź dla trzech godzin. Jeśli występują dla trzech godzin to sprawdzono, że co najwyżej dwa braki występują po sobie, więc zdecydowano się je zastąpić średnią z 2 godzin poprzedzających i następujących po braku.

In [34]:
(data_merged['PM2.5'].isna().groupby(data_merged['Kod stacji']).mean() * 100).sort_values(ascending=False)
Out[34]:
Kod stacji
DsDusznikMOB    0.0
PkPrzemGrunw    0.0
MzWarKondrat    0.0
MzWarTolstoj    0.0
MzWarWokalna    0.0
               ... 
MzOtwoBrzozo    0.0
MzPiasPulask    0.0
MzPlocMiReja    0.0
MzRadTochter    0.0
ZpSzczPilsud    0.0
Name: PM2.5, Length: 62, dtype: float64

Sprawdzenie,czy wszystkie dane nie zawierają już brakujących wartości PM2.5.

In [35]:
plt.figure(figsize=(12, 6))

for kod_stacji in data_merged['Kod stacji'].unique():
    dane_stacji = data_merged[data_merged['Kod stacji'] == kod_stacji]
    plt.plot(dane_stacji['Date'], dane_stacji['PM2.5'])

plt.title('PM2.5 dla każdej stacji')
plt.xlabel('Data i godzina')
plt.ylabel('PM2.5')
plt.xticks(rotation=45)
plt.show()

TWORZENIE DODATKOWYCH CECH Z DATATIME¶

In [36]:
data_merged['Date'] = pd.to_datetime(data_merged['Date'])
data_merged['day'] = data_merged['Date'].dt.dayofweek
data_merged['month'] = data_merged['Date'].dt.month
data_merged['weekend'] = data_merged['day'].isin([5, 6]).astype(int)
data_merged['year'] = data_merged['Date'].dt.year

def get_season(month):
    if month in [3, 4, 5]:
        return 0
    elif month in [6, 7, 8]:
        return 1
    elif month in [9, 10, 11]:
        return 2
    else:
        return 3

data_merged['season'] = data_merged['month'].apply(get_season)

Zdecydowano się dodać nowe cechy takie jak: dzień, miesiąc, rok oraz informację czy jest weekend czy dzień roboczy.

CECHY W OKNACH CZASOWYCH¶

In [37]:
unique_station_ids = data_merged['Kod stacji'].unique()
random_station_ids = random.sample(list(unique_station_ids), 1)
print(random_station_ids)
['LbLubObywate']
In [38]:
stacja = data_merged[(data_merged['month'] == 1) & (data_merged['year'] == 2019)]
stacja = stacja[stacja['Kod stacji'] == 'WmOlsPuszkin']
stacja = stacja.set_index("Date")

statistics = ["mean", "std", "max"]
windows = ['1h', '3h', '6h', '12h','1D', '3D']

for window in windows:
    rolling_df = stacja['PM2.5'].rolling(window = window, closed='left').agg(statistics)
    fig, axes = plt.subplots(len(statistics), 1, figsize = [20, 16])
    fig.suptitle(f'Okno = {window}')

    for idx, stat in enumerate(statistics):
        sns.lineplot(data=stacja['PM2.5'], ax = axes[idx], color='red', label='PM2.5')
        sns.lineplot(data=rolling_df[stat], ax=axes[idx], color='blue', label=f'{stat} PM2.5')
        axes[idx].set_ylabel(f'{stat.capitalize()} PM2.5')

plt.show()
In [39]:
# stacja = data_merged[(data_merged['month'] == 1) & (data_merged['year'] == 2019)]
data_stacja = data_merged[data_merged['Kod stacji'] == 'WmOlsPuszkin']
data_stacja = data_stacja.set_index("Date")

statistics = ["mean", "std", "max"]
windows = ['7D', '30D']

for window in windows:
    rolling_df = data_stacja['PM2.5'].rolling(window = window, closed='left').agg(statistics)
    fig, axes = plt.subplots(len(statistics), 1, figsize = [20, 16])
    fig.suptitle(f'Okno = {window}')

    for idx, stat in enumerate(statistics):
        sns.lineplot(data=data_stacja['PM2.5'], ax = axes[idx], color='red', label='PM2.5')
        sns.lineplot(data=rolling_df[stat], ax=axes[idx], color='blue', label=f'{stat} PM2.5')
        axes[idx].set_ylabel(f'{stat.capitalize()} PM2.5')

plt.show()

Od razu widać, że jednogodzinne okno nie niesie żadnej wartości, analogicznie dla 3godzinnego oraz 6godzinnego okna. Średnia dla 12 godzin ma już jakikolwiek sens. Zdecydowano się wybrać również okno ze średnią 1dniową oraz maximum 1dniowe. Wybrano również maxmium 30dniowe.

In [40]:
data_merged.set_index('Date', inplace=True)

def calculate_rolling_stats(group):
    group['rolling_mean_12_h'] = group['PM2.5'].rolling(window='12h', closed='left').mean()
    group['rolling_mean_1D'] = group['PM2.5'].rolling(window='1D', closed='left').mean()
    group['rolling_max_1D'] = group['PM2.5'].rolling(window='1D', closed='left').max()
    group['rolling_max_30D'] = group['PM2.5'].rolling(window='30D', closed='left').max()
    return group

data_with_rolling_stats = data_merged.groupby('Kod stacji').apply(calculate_rolling_stats)
data_with_rolling_stats.reset_index(inplace=True)
C:\Users\User\AppData\Local\Temp\ipykernel_23260\2148299505.py:10: FutureWarning: Not prepending group keys to the result index of transform-like apply. In the future, the group keys will be included in the index, regardless of whether the applied function returns a like-indexed object.
To preserve the previous behavior, use

	>>> .groupby(..., group_keys=False)

To adopt the future behavior and silence this warning, use 

	>>> .groupby(..., group_keys=True)
  data_with_rolling_stats = data_merged.groupby('Kod stacji').apply(calculate_rolling_stats)

Dodano do zbioru danych cechy tj. średnia 12 godzinna, średnia 1 dniowa, maximum 1 dniowe oraz maximum 30 dniowe.

CECHY OPARTE NA OPÓŹNIENIACH¶

In [41]:
unique_station_ids = data_merged['Kod stacji'].unique()
random_station_ids = random.sample(list(unique_station_ids), 5)
print(random_station_ids)
['KpBydPlPozna', 'WmGoldUzdrowMOB', 'MzWarChrosci', 'PkMielBierna', 'MzZyraRoosev']
In [42]:
data_jan_2020 = data_merged[(data_merged['year'] == 2019) & (data_merged['month'] == 1)]

for i in random_station_ids:
    pm25 = data_jan_2020[data_jan_2020['Kod stacji'] == i]['PM2.5']
    plot_pacf(pm25, lags=24)
    plt.show()

Dla większości przypadów najsilniejszą korelację osiągają opóźnienia 1godzinne, 2godzinne, 3godzinne, 20godzinne oraz 24 godzinne. Dlatego też zdecydowano się takie lagi dodać.

In [43]:
def calculate_rolling_stats(group):
    lag1 = group['PM2.5'].shift(1)
    group['lag1_h'] = lag1

    lag2 = group['PM2.5'].shift(2)
    group['lag2_h'] = lag2

    lag3 = group['PM2.5'].shift(3)
    group['lag3_h'] = lag3

    lag20 = group['PM2.5'].shift(20)
    group['lag20_h'] = lag20

    lag24 = group['PM2.5'].shift(24)
    group['lag24_h'] = lag24
    return group

data_with_lags = data_with_rolling_stats.groupby('Kod stacji').apply(calculate_rolling_stats)
data_with_lags.reset_index(inplace=True)
data_with_lags = data_with_lags.drop(['index'], axis = 1)
C:\Users\User\AppData\Local\Temp\ipykernel_23260\2090405549.py:18: FutureWarning: Not prepending group keys to the result index of transform-like apply. In the future, the group keys will be included in the index, regardless of whether the applied function returns a like-indexed object.
To preserve the previous behavior, use

	>>> .groupby(..., group_keys=False)

To adopt the future behavior and silence this warning, use 

	>>> .groupby(..., group_keys=True)
  data_with_lags = data_with_rolling_stats.groupby('Kod stacji').apply(calculate_rolling_stats)
In [44]:
data_with_lags.head(5)
Out[44]:
Date Kod stacji PM2.5 Nazwa stacji Typ obszaru Województwo Miejscowość lat lon srednia_wojewodzka_godz ... season rolling_mean_12_h rolling_mean_1D rolling_max_1D rolling_max_30D lag1_h lag2_h lag3_h lag20_h lag24_h
0 2019-01-01 01:00:00 DsDusznikMOB 33.40530 Duszniki-Zdrój miejski DOLNOŚLĄSKIE Duszniki-Zdrój 50.402645 16.393319 82.543420 ... 3 NaN NaN NaN NaN NaN NaN NaN NaN NaN
1 2019-01-01 02:00:00 DsDusznikMOB 13.80280 Duszniki-Zdrój miejski DOLNOŚLĄSKIE Duszniki-Zdrój 50.402645 16.393319 54.379220 ... 3 33.405300 33.405300 33.4053 33.4053 33.40530 NaN NaN NaN NaN
2 2019-01-01 03:00:00 DsDusznikMOB 9.94056 Duszniki-Zdrój miejski DOLNOŚLĄSKIE Duszniki-Zdrój 50.402645 16.393319 43.740992 ... 3 23.604050 23.604050 33.4053 33.4053 13.80280 33.40530 NaN NaN NaN
3 2019-01-01 04:00:00 DsDusznikMOB 6.75889 Duszniki-Zdrój miejski DOLNOŚLĄSKIE Duszniki-Zdrój 50.402645 16.393319 34.282634 ... 3 19.049553 19.049553 33.4053 33.4053 9.94056 13.80280 33.4053 NaN NaN
4 2019-01-01 05:00:00 DsDusznikMOB 7.88722 Duszniki-Zdrój miejski DOLNOŚLĄSKIE Duszniki-Zdrój 50.402645 16.393319 22.843292 ... 3 15.976887 15.976887 33.4053 33.4053 6.75889 9.94056 13.8028 NaN NaN

5 rows × 24 columns

In [45]:
data_with_lags.set_index('Date', inplace = True)
data_with_lags.sort_index(inplace = True)

IMPORTY BIBLIOTEK¶

In [46]:
from sklearn.preprocessing import LabelEncoder
from sklearn.model_selection import train_test_split, GridSearchCV
import xgboost as xgb
from sklearn.metrics import mean_squared_error, r2_score, mean_squared_error
from sklearn.model_selection import TimeSeriesSplit

INICJALIZACJA I TRENOWANIE MODELU¶

In [47]:
X = data_with_lags.drop(['PM2.5'], axis=1)
y = data_with_lags['PM2.5']
In [48]:
columns_to_encode = [0, 1, 2, 3, 4]

for c in columns_to_encode:
    labelencoder = LabelEncoder()
    X.iloc[:, c] = labelencoder.fit_transform(X.iloc[:, c])
C:\Users\User\AppData\Local\Temp\ipykernel_23260\3553582975.py:5: DeprecationWarning: In a future version, `df.iloc[:, i] = newvals` will attempt to set the values inplace instead of always setting a new array. To retain the old behavior, use either `df[df.columns[i]] = newvals` or, if columns are non-unique, `df.isetitem(i, newvals)`
  X.iloc[:, c] = labelencoder.fit_transform(X.iloc[:, c])
C:\Users\User\AppData\Local\Temp\ipykernel_23260\3553582975.py:5: DeprecationWarning: In a future version, `df.iloc[:, i] = newvals` will attempt to set the values inplace instead of always setting a new array. To retain the old behavior, use either `df[df.columns[i]] = newvals` or, if columns are non-unique, `df.isetitem(i, newvals)`
  X.iloc[:, c] = labelencoder.fit_transform(X.iloc[:, c])
C:\Users\User\AppData\Local\Temp\ipykernel_23260\3553582975.py:5: DeprecationWarning: In a future version, `df.iloc[:, i] = newvals` will attempt to set the values inplace instead of always setting a new array. To retain the old behavior, use either `df[df.columns[i]] = newvals` or, if columns are non-unique, `df.isetitem(i, newvals)`
  X.iloc[:, c] = labelencoder.fit_transform(X.iloc[:, c])
C:\Users\User\AppData\Local\Temp\ipykernel_23260\3553582975.py:5: DeprecationWarning: In a future version, `df.iloc[:, i] = newvals` will attempt to set the values inplace instead of always setting a new array. To retain the old behavior, use either `df[df.columns[i]] = newvals` or, if columns are non-unique, `df.isetitem(i, newvals)`
  X.iloc[:, c] = labelencoder.fit_transform(X.iloc[:, c])
C:\Users\User\AppData\Local\Temp\ipykernel_23260\3553582975.py:5: DeprecationWarning: In a future version, `df.iloc[:, i] = newvals` will attempt to set the values inplace instead of always setting a new array. To retain the old behavior, use either `df[df.columns[i]] = newvals` or, if columns are non-unique, `df.isetitem(i, newvals)`
  X.iloc[:, c] = labelencoder.fit_transform(X.iloc[:, c])

Zenkodowano kolumny kategoryczne, aby móc użyć je później w modelu.

In [49]:
X.head(5)
Out[49]:
Kod stacji Nazwa stacji Typ obszaru Województwo Miejscowość lat lon srednia_wojewodzka_godz day month ... season rolling_mean_12_h rolling_mean_1D rolling_max_1D rolling_max_30D lag1_h lag2_h lag3_h lag20_h lag24_h
Date
2019-01-01 01:00:00 0 4 0 0 4 50.402645 16.393319 82.543420 1 1 ... 3 NaN NaN NaN NaN NaN NaN NaN NaN NaN
2019-01-01 01:00:00 5 3 0 1 3 53.121764 17.987906 46.479233 1 1 ... 3 NaN NaN NaN NaN NaN NaN NaN NaN NaN
2019-01-01 01:00:00 58 11 0 10 12 51.747950 18.049063 75.200000 1 1 ... 3 NaN NaN NaN NaN NaN NaN NaN NaN NaN
2019-01-01 01:00:00 43 7 1 7 8 50.192100 23.361425 40.202483 1 1 ... 3 NaN NaN NaN NaN NaN NaN NaN NaN NaN
2019-01-01 01:00:00 24 22 0 4 21 52.178766 21.560968 45.312354 1 1 ... 3 NaN NaN NaN NaN NaN NaN NaN NaN NaN

5 rows × 22 columns

In [50]:
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25, shuffle = False)
In [51]:
fig, axes = plt.subplots(1,1, figsize = (20,10))
plt.plot(y_train, color = 'red')
plt.plot(y_test, color = 'blue')
Out[51]:
[<matplotlib.lines.Line2D at 0x14fafb94b90>]

Podzielono zbiór danych na treningowy i testowy w stosunku 75% do 25%. Zdecydowano się na wykonanie predykcji na danych ogólnych (nie dla każdej stacji z osobna) ze względu, że jest to bardziej optymalne obliczeniowo i mamy jeden model, który możemy zastosować dla całości danych, a nie musimy zależnie od stacji używać innego modelu.

In [52]:
model = xgb.XGBRegressor()

model.fit(X_train, y_train)
y_pred = model.predict(X_test)

rmse = mean_squared_error(y_test, y_pred, squared=False)
print("RMSE:", rmse)
RMSE: 6.247570667847238

Błąd średniokwadratowy dla modelu wynosi około 6.25% przy domyślnych parametrach.

PLUSY I MINUSY WYKORZYSTANIA XGB¶

PLUSY:

+zawiera techniki regularyzacji, które pomagają zapobiec nadmiernemu dopasownaniu,

+wysoka dokładność,

+dobrze radzi sobie z dużymi datasetami,

+nie ma potrzeby normalizowania danych, radzi sobie z brakującymi wartościami

MINUSY:

-wymaga dostosowania hiperparametróww, przekonano się, że jest to bardzo czasochłonne,

-pomimo technik regularyzacji algorytm nadal może przeuczyć dane treningowe,

-interpretacja i zrozumienie w jaki sposób algorytm osiąga przewidywania jest trudna,

-trudne do dostrojenia ze względu na dużą ilość hiperparametrów.

TUNING HIPERPARAMETRÓW¶

In [53]:
param_grid = {
    'max_depth': [3, 5, 7],
    'learning_rate': [0.1, 0.01, 0.001],
    'n_estimators': [100, 200, 300],
}

grid_search = GridSearchCV(estimator=xgb.XGBRegressor(), param_grid=param_grid, cv=3, scoring='neg_mean_squared_error')
grid_search.fit(X_train, y_train)

print("Najlepsze parametry:", grid_search.best_params_)

best_model = grid_search.best_estimator_
y_pred = best_model.predict(X_test)
Najlepsze parametry: {'learning_rate': 0.1, 'max_depth': 5, 'n_estimators': 100}

Po wykonaniu tuningu parametrów najlepszymi okazują się być:

-'learning_rate': 0.1,

-'max_depth': 5,

-'n_estimators': 100

Próbowano wykonywać kolejne tuningi z większą liczbą parametrów, lecz procesy te były zbyt czasochłonne.

CZYM SĄ HIPERPARAMETRY¶

Hiperparametry to zewnętrzne zmienne konfiguracyjne, których Data Scientiści używają do zarządzania szkoleniem modeli uczenia maszynowego. Są ustawiane ręcznie przed szkoleniem modelu. Hiperparametry określają kluczowe funkcje, takie jak architektura modelu, szybkość uczenia się i złożoność modelu. Mają znaczący wpływ na dokładność i skuteczność modelu.

Strojenie hiperparametrów to proces wybierania optymalnych wartości hiperparametrów modelu uczenia maszynowego. Celem dostrajania hiperparametrów jest znalezienie wartości, które prowadzą do najlepszej wydajności danego zadania.

In [54]:
best_params = {'learning_rate': 0.1, 'max_depth': 5, 'n_estimators': 100}
final_model = xgb.XGBRegressor(**best_params)

final_model.fit(X_train, y_train)

y_pred = final_model.predict(X_test)

rmse = mean_squared_error(y_test, y_pred, squared=False)
print("RMSE:", rmse)
RMSE: 6.198891964836731

Po dostrojeniu hiperparametrów błąd średniokwadratowy zmalał w porównaniu do modelu z wartościami domyślnymi.

In [55]:
testing_data = pd.DataFrame(data=X_test, columns=X.columns)
testing_data['PM2.5'] = y_test  
testing_data['pm25_przewidywane'] = y_pred
testing_data.head()
Out[55]:
Kod stacji Nazwa stacji Typ obszaru Województwo Miejscowość lat lon srednia_wojewodzka_godz day month ... rolling_mean_1D rolling_max_1D rolling_max_30D lag1_h lag2_h lag3_h lag20_h lag24_h PM2.5 pm25_przewidywane
Date
2019-10-01 19:00:00 33 48 0 4 41 52.290864 21.042458 7.482966 1 10 ... 4.482526 6.322227 51.180685 5.011755 4.507727 4.507727 4.709338 1.987589 5.213366 5.518918
2019-10-01 19:00:00 26 28 1 4 27 52.191728 20.837489 7.482966 1 10 ... 5.584167 7.210000 260.490000 6.190000 5.400000 5.510000 5.060000 4.610000 7.100000 6.485516
2019-10-01 19:00:00 58 11 0 10 12 51.747950 18.049063 NaN 1 10 ... 6.613928 11.300000 65.493700 4.094110 4.182160 2.128030 8.866600 5.002220 7.516842 7.554872
2019-10-01 19:00:00 19 16 0 5 16 50.010575 19.949189 9.795942 1 10 ... 5.686784 15.587000 52.089600 6.675240 3.294000 3.873900 4.179670 8.146100 5.802210 7.633582
2019-10-01 19:00:00 9 20 0 2 19 51.259431 22.569133 8.867390 1 10 ... 4.987500 8.800000 64.600000 8.800000 7.700000 5.600000 4.400000 2.900000 11.900000 9.025276

5 rows × 24 columns

In [56]:
testing_data_stacja = testing_data[testing_data['Kod stacji'] == 33]
r2 = r2_score(y_test, y_pred)

print(f'R-squared: {r2}')
plt.figure(figsize=(12, 8))
plt.plot(testing_data_stacja['PM2.5'], label='Rzeczywiste PM2.5')
plt.plot(testing_data_stacja['pm25_przewidywane'], label='Przewidywane PM2.5')

plt.legend()
plt.title('Rzeczywiste vs przewidywane wartości PM2.5')
plt.xlabel('Index')
plt.ylabel("Stężenie PM2.5 (µg/m³)")
plt.show()
R-squared: 0.8744843687608252
In [57]:
# param_grid = {
#     'max_depth': [3, 5, 7],
#     'learning_rate': [0.05, 0.1, 0.2],
#     'n_estimators': [100, 200, 300],
#     'colsample_bytree': [0.6, 0.8, 1.0],
#     'subsample': [0.6, 0.8, 1.0],
#     'gamma': [0, 0.1, 0.2],
#     'reg_alpha': [0, 0.1, 0.2],
#     'reg_lambda': [0, 0.1, 0.2],
#     'min_child_weight': [1, 3, 5],
#     'max_delta_step': [0, 1, 2]
# }

# grid_search = GridSearchCV(estimator=xgb.XGBRegressor(), param_grid=param_grid, cv=3, scoring='neg_mean_squared_error')

# grid_search.fit(X_train, y_train)

# print("Najlepsze parametry:", grid_search.best_params_)

# best_model = grid_search.best_estimator_
# y_pred = best_model.predict(X_test)

# rmse = mean_squared_error(y_test, y_pred, squared=False)
# print("RMSE:", rmse)

BACKTESTING¶

In [58]:
tscv = TimeSeriesSplit(n_splits=5)

Walidacja Walk Forward (WFV) to technika walidacji krzyżowej szeregów czasowych stosowana do oceny wydajności modeli predykcyjnych. Walk forward polega na podziale zbioru danych na sekwencje ciągłych okien czasowych. Jest to szczególnie przydatne w przypadku danych uporządkowanych w czasie, gdzie sekwencja danych ma znaczenie. WFV zaprojektowano tak, aby był bardziej realistyczny w ocenie, jak dobrze model będzie generalizował na przyszłe, niewidoczne dane.

DZIAŁANIE:

Początkowy okres szkolenia i testu.

Trenowanie modelu.

Model testowy.

Przesunięcie okna.

Powtórzenie.

In [59]:
rmse_scores = []

for train_index, test_index in tscv.split(X_train):
    X_train_fold, X_test_fold = X_train.iloc[train_index], X_train.iloc[test_index]
    y_train_fold, y_test_fold = y_train.iloc[train_index], y_train.iloc[test_index]
    print(f"Długości danych treningowych: {len(y_train_fold)}")
    model.fit(X_train_fold, y_train_fold)
    y_pred = model.predict(X_test_fold)
    
    rmse = mean_squared_error(y_test_fold, y_pred, squared=False)
    rmse_scores.append(rmse)

mean_rmse = sum(rmse_scores) / len(rmse_scores)
print("Mean RMSE:", mean_rmse)
Długości danych treningowych: 67890
Długości danych treningowych: 135780
Długości danych treningowych: 203670
Długości danych treningowych: 271560
Długości danych treningowych: 339450
Mean RMSE: 3.8984356076768236

Czym rózni się backtesting od testowania wydajności modelu dla danych niebędących szeregami czasowymi?¶

Backtesting polega na ocenie wydajności modelu na danych historycznych, które nie zostały użyte do trenowania, symulując prawdziwe warunki na rynku lub środowisku. Celem backtestingu jest sprawdzenie, jak dobrze model radzi sobie ze zrozumieniem i przewidywaniem rzeczywistych zdarzeń na podstawie historycznych danych.

Testowanie wydajności modelu dla danych niebędących szeregami czasowymi odbywa się na zbiorze danych, który jest jednorodny w czasie, ale nie ma konkretnego związku z określonymi punktami w czasie. Celem testowania wydajności modelu dla danych niebędących szeregami czasowymi jest ocena ogólnej zdolności modelu do generalizacji na nowych, nieznanych danych, które nie mają struktury czasowej.

Podsumowując, backtesting ma specyficzny charakter związany z oceną modeli predykcyjnych w kontekście historycznych danych czasowych, podczas gdy testowanie wydajności modelu dla danych niebędących szeregami czasowymi koncentruje się na ocenie zdolności modelu do generalizacji na różne rodzaje danych bez struktury czasowej.

Analiza istotności cech¶

In [60]:
feature_importance = final_model.feature_importances_

feature_importance_df = pd.DataFrame({'Feature': X_train.columns, 'Importance': feature_importance})

feature_importance_df = feature_importance_df.sort_values(by='Importance', ascending=False)

plt.figure(figsize=(10, 6))
plt.barh(feature_importance_df['Feature'], feature_importance_df['Importance'])
plt.xlabel('Importance')
plt.ylabel('Feature')
plt.title('Feature Importance')
plt.show()
In [61]:
feature_importance_df
Out[61]:
Feature Importance
17 lag1_h 0.836431
7 srednia_wojewodzka_godz 0.048263
15 rolling_max_1D 0.030158
16 rolling_max_30D 0.009065
1 Nazwa stacji 0.008964
19 lag3_h 0.007767
8 day 0.007317
18 lag2_h 0.005483
4 Miejscowość 0.005044
6 lon 0.004881
5 lat 0.004796
21 lag24_h 0.004277
20 lag20_h 0.004110
12 season 0.003913
13 rolling_mean_12_h 0.003774
14 rolling_mean_1D 0.003760
2 Typ obszaru 0.003722
0 Kod stacji 0.003030
9 month 0.002894
3 Województwo 0.002350
10 weekend 0.000000
11 year 0.000000

Można zaobserować (tak jak można było się spodziewać), że największą istotność ma kolumna z wartościami jednogodzinnego opóźnienia. Kolejnymi czynnikami, które mają jakiekolwiek znaczenie są cechy stworzone przy okazji podziału na okna czasowe oraz w trakcie tworzenia opóźnień.